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We study the attractive fermionic Hubbard model on a honeycomb lattice using determinantal 
quantum Monte Carlo simulations. By increasing the interaction strength U (relative to the hopping 
parameter t) at half-filling and zero temperature, the system undergoes a quantum phase transition 
at 5.0 < Uc/t < 5.1 from a semi-metal to a phase displaying simultaneously superfluid behavior and 
density order. Doping away from half-filling, and increasing the interaction strength at finite but low 
temperature T, the system always appears to be a superfluid exhibiting a crossover between a BCS 
and a molecular regime. These different regimes are analyzed by studying the spectral function. 
The formation of pairs and the emergence of phase coherence throughout the sample are studied as 
U is increased and T is lowered. 
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The recent discovery of graphene layers, i.e. single- 
atom thick layers of carbon atoms arranged in a planar 
honeycomb structure has attracted considerable atten- 
tion due to its interest in fundamental physics as well 
as for potential applications. The energy band spectrum 
shows "conical points" where the valence and conduc- 
tion bands are connected, and the Fermi energy at half- 
filling is located precisely at these points as only half 
of the available states arc filled. Around these points, 
the energy varies proportionally to the modulus of the 
wave- vector and the excitations (holes or particles) of 
the system are equivalent to ultra-rclativistic (massless) 
Dirac fermions since their dispersion relation is lineari^ 
Graphene sheets then allows for table-top experiments on 
two-dimensional field theories with quantum anomalies, 
allowing us to explore the Klein paradox;^ the anomalous 
quantum Hall effect induced by Berry phases^ and its 
corresponding modified Landau levelsi^ 

When the fermions are interacting, the peculiar na- 
ture of the Fermi surface (i.e. reduced to a finite number 
of Dirac points) leads to special physics at and around 
half-filling. In a square lattice, the nesting of the Fermi 
surface generally leads to ordered phases even for arbi- 
trarily small interaction strengths. On the contrary, in 
the honeycomb lattice and with repulsive interactions, 
Paiva et al. have found^ a quantum phase transition 
(QPT) at half-filling between a metallic and an ordered 
phase when the interaction strength is increased. How- 
ever, since graphene is a weakly-interacting system, this 
QPT is not accessible experimentally. 

In a recent work, some of us have analyzed the possibil- 
ity of reproducing graphene physics and of extending it 
to the interacting regime by creating a two-dimensional 
honeycomb optical lattice and loading ultracold spin- 1/2 
fermionic atoms, such as ^Li, into it4 The key advantage 



is that the relevant experimental parameters {e.g. config- 
uration and strength of the optical potential, inter-atomic 
interaction strength tuned via Feshbach resonance) can 
be accurately controlled while getting rid of the inherent 
complexity of a solid. Following this idea, we use ex- 
act Quantum Monte Carlo (QMC) simulations to study 
interacting ultracold fermions loaded into a honeycomb 
optical lattice in the absence of any external confinement. 
We will focus on the case of attractive interactions as it is 
accessible with these numerical techniques and free from 
the sign problem at and away from half-filling. 

In the continuum at zero temperature, as the inter- 
acting fermionic gas is driven from the weak to the 
strong attractive coupling limit, there is a crossover from 
a BCS regime of weakly-bound delocalized pairs to a 
Bose-Einstein condensate (BEC) of tightly-bound pairs 
(later called molecules for simplicity) ^° i ^ ^ At finite 
but sufficiently low temperature, a similar BCS-molecule 
crossover is observed except that, the system being two- 
dimensional, there is only quasi-long-range order and, 
consequently, no true condensate but only a superfluid. 
In this paper, we will study interacting particles on a lat- 
tice, represented by a simple fermionic Hubbard modelJ^ 
Nonetheless, some aspects of the continuum limit, such 
as the BCS-BEC crossover, are expected to be repro- 
duced in the discrete model. Zhao and Paramekanti have 
explored the attractive fermionic Hubbard model on a 
honeycomb lattice using mean field theoryi^ and they 
found a QPT between a semi-metal and a superfluid at 
half-filling. Away from half-filling, they recovered the 
crossover already observed in the continuum limit. Re- 
cently, Su et al. used QMC methods to study the BCS- 
BEC crossover on the honeycomb lattice away from half- 
filling and concluded that it was similar to the one ob- 
tained for the square latticCf^ In the present work, we 
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use QMC simulations and large system sizes to study 
the pair formation at half-filling and accurately deter- 
mine the critical value of the coupling strength at which 
pairs form. We then study pairing away from half-filling 
by analyzing several quantities, including spectral func- 
tions. 

The paper is organized as follows. In section I, we in- 
troduce the model, notations and the quantities we use to 
characterize the different phases. In section II, we show 
that our system at half-filling can be related to the repul- 
sive Hubbard model^ and then present complementary 
results for this case, including the QPT point the system 
crosses to go from a semi-metallic disordered phase to 
an ordered one displaying both superfluid behavior and 
density wave order. The location of this QPT point has 
been accurately determined compared to previous works, 
and the nature of the weakly-interacting phase before the 
transition is addressed by analyzing the behavior of the 
spectral function as the interaction strength is varied. 
Finally, in section III we study the system doped away 
from half-filling. The system is clearly shown to exhibit 
superfluid behavior while the density wave order present 
at half-filling has been destroyed. We conclude our study 
by analyzing the formation of pairs and the emergence of 
global phase coherence as a function of temperature and 
interaction strength. 



are massless Dirac fermions that obey the 2D Wcyl-Dirac 
equationii^ 

The FAHM ([I]) on a bipartite lattice is particle-hole 
symmetrioii and thus adopts the same phases for densi- 
ties p and 2 — p. It is then sufRcicnt to study the system 
for densities p > I. This model can also be mapped onto 
the fermionic repulsive Hubbard model (FRHM)iii^ by 
performing a particle-hole transformation on only one of 
the species. Consequently, the physics of the FAHM at 
densities (pfiPi) is equivalent to that of the FRHM at 
densities (1 — pj, pj^) or 1 — p{), but with a non-zero 
Zeeman-like term, — /iX^i ~ ^ii)- Therefore, the two 
models are identical at half- filling [p, = 0). We will use 
this equivalence in section|lT]where we concentrate on the 
half- filled case. 




I. THE FERMIONIC HUBBARD MODEL 

The physics of a system of A'f spin- 1/2 fermions, with 
attractive two-body interactions and equal spin popula- 
tions, filling up a lattice made of sites is encapsulated 
in a simple tight-binding model, namely the fermionic at- 
tractive Hubbard model (FAHM) , whose grand-canonical 
Hamiltonian operator readsJ^ 



(1) 



Here denotes pairs of nearest-neighbors sites on 

the lattice, a =tji a-re the two possible spin states of 
the fermions, /^^ and fia- are the creation and annihila- 
tion operators of a fermion with spin state a at site i, 
"■io- = fiafiiy is the corresponding number operator, t is 
the hopping amplitude between nearest-neighbors sites, 
{/ > is the strength of the attractive interaction be- 
tween fermions with opposite spin states and p is the 
chemical potential whose value fixes the average total 
fermionic density p. With the present form of the in- 
teraction term, the system is half- filled, i.e. there is on 
average one fermion per site [p = Nf/N = I), when 
/i = 0. In the non-interacting limit [7 = 0, this system is 
known to behave like a semi-metal with vanishing density 
of states at the Fermi level and its elementary excitations 



FIG. 1: Finite honeycomb lattice of linear dimension L = 3. 
The total number of sites is = 2L^ = f8. 



To calculate the equilibrium properties of this model 
at finite but low temperatures T, we used the stan- 
dard determinant quantum Monte Carlo algorithm 
(DQMC)3i^'^i2ii^ The cases under our consideration 
(namely attractive interactions and equal densities of 
spin-up and spin-down fermions) are free of the sign 
problem^! that used to plague numerical simulations of 
fermionic systems. This will allow us to reach the low 
temperatures needed to study pairing and superfluidity. 
In the following, the reciprocal of the thermal energy 
(also called the inverse temperature) is denoted as usual 
by /3 = l/ksT, where fcs is the Boltzmann constant. 

In the DQMC simulations, we have used the honey- 
comb lattice depicted in Fig. [l]with periodic boundary 
conditions. The primitive vectors ai and a.2 delineate 
a diamond-shaped primitive cell of the Bravais lattice 
which contains two nonequivalent sites (a and b) sepa- 
rated by AB = (ai +a2)/3 and each producing upon tiling 
a hexagonal sublattice. A finite honeycomb lattice of side 
L then contains N ~ 2L^ sites. In the non-interacting 
case, the energy levels are given hy^^ 



e±(fci,fc2) =±i 1 



where fci,/c2 € {0, 1, • • • , i — 1}. When L is a multi- 
ple of three, there always exist pairs (fci,fc2) such that 
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FIG. 2: (Color online) Total average density p vs chemical 
potential fi for U/t — (top) and U/t — 1 (bottom) at /3f = 16 
and different lattice sizes L. The top figure is obtained by 
analytical calculation at (7 = 0. The bottom figure is obtained 
fi'om numerical data generated by DQMC. For sizes that are 
not multiples of three, there is no state at half-filling and a 
small gap appears for small system sizes. There is no such gap 
when L is a multiple of three. For sizes that are multiples of 
three, plateaus appear away from half-filling. These plateaus 
are also finite-size effects and they disappear when L —> oo. 
The dotted line in the top figure is obtained by an exact 
evaluation of the derivative dp/dp\^=o in the non-interacting 
limit when L — > oo. The two figures show that the "magic 
number 3" effect is present even when the interaction strength 
U is comparable to the hopping parameter t. 



e±{ki,k2) — 0, i.e. there are four states (two per spin 
state) located exactly at the Fermi level and only two of 
these states will be occupied if p = 1 . This does not hap- 
pen when L is not a multiple of three. As a consequence, 
on small finite-size systems, a small gap of order 1/L ap- 
pears around half-filling when L is not a multiple of three 
(see Fig. To avoid confusion between this gap, which 
is a finite-size effect, and Mott gaps generated by interac- 
tions that are expected to appear in ordered phases, we 



used (especially at half-filling) sizes L that are multiples 
of three. This limits strongly the sizes that can be stud- 
ied, fn the most favorable cases, we went up to L = 15, 
that is iV = 450 sites. 

In the strong coupling regime (U S> i), we expect 
the system to form pairs (hereafter called molecules) of 
fcrmions with opposite spins on the same site. These 
pairs can show two different ordering phenomena: estab- 
lishment of a phase coherence order or of a solid (crystal- 
type) order. A solid of pairs would exhibit a density wave 
typical of a crystal and would reveal itself through spatial 
oscillations in the density-density correlation function. 



(2) 



where rii ~ ^i<^ '^^ ^^'^ total number of fermions on site 
i and where (•) denotes the quantum statistical average at 
temperature T . At half-filling and zero temperature, we 
expect to observe a phase where alternate sites are empty 
and where only the A or the B sub-lattice is occupied. 
Such a density wave is signaled by a structure factor 5'dw 
diverging linearly with the total number of sites TV of the 
system, where 



S'dv 



(3) 



with the site index i being even on A sites and odd on B 
sites. 

In a Bose condensed phase, the phase coherence be- 
tween pairs is signaled by long-range order (or quasi-long- 
range order for a superfluid at finite temperature) in the 
pair Green's function. 



(4) 



where Aj — fl-^fly creates a pair on site i. In a way simi- 
lar to the density correlations, we define a pair structure 
factor P^r^ 



(5) 



This pair structure factor diverges linearly with N when 
long-range order is achieved. Finally, in the absence of 
any order, the system is expected to be a semi-metal at 
half-filling due to the peculiar nature of the Fermi surface 
(no gap but a vanishing density of states at the Fermi 
level). To distinguish between metallic, semi- metallic 
or gapped (solid or superfluid) states, we calculate the 
spectral function A{ijj) which essentially reflects the one- 
particle density of states. To obtain this quantity, we first 
calculate the (imaginary) time-displaced on-site Green's 
function G{t) = T,t{ f ti^) fKo)) /N and then extract 
A{lu) by inverting the following Laplace transform 



G(r) 



duj- 



using an analytic continuation methodi^ 
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II. HONEYCOMB LATTICE AT HALF-FILLING 

At half-filling, the system can be mapped onto the 
FRHMi^i^i2^i2i Defining a hole creation operator h]^ 
for the down spin through, 

(-ir/»l^ = /a, (6) 

the kinetic term is left unchanged in the spin-down holes 
representation. The number operator is accordingly 
transformed into 1 — n^, where n^j^ — ^1x^4 is the num- 
ber operator for holes, and, up to a redefinition of the 
chemical potential fi, the sign of the interaction term is 
reversed. The FRHM has SU(2) spin-rotation symmetry 
at half-filling, which translates into the SU(2) pseudo- 
spin symmetry of FAHMi^ Hence the spin-spin correla- 
tions are the same along the three coordinate axes. 

{ata]) = {ayay)^{ata^), (7) 

where x and y are the in-plane axes and z the axis or- 
thogonal to the lattice plane. More specifically: 

< = fl^Ki + h^af^^, (8) 

At large interaction, the FRHM is known to be equiva- 
lent to a Heisenberg model and it develops a long-range 
anti-ferromagnetic order on the honeycomb lattice at zero 
temperatureii The correlation functions ([7]) then show 
oscillations from site to site. Translated into the attrac- 
tive model language, these functions bccom e^^i^^ 

= {n,nj~n,-nj-\), (9) 
«'^J + '^H> = 2(-iy+MAlA,+A,At). (10) 

The spin anti-fcrromagnctic correlations along the z-axis 
in the FRHM are then reproduced in the density-density 
correlations Dij of the FAHM, which develops a den- 
sity wave with alternating occupied and empty sites. 
The spin correlations in the xy lattice plane translate 
into long-range order for the Green's function G?- and 
phase coherence of a Bose-Einstein condensate. The anti- 
ferromagnetic phase of the FRHM is thus mapped onto 
a peculiar phase for the FAHM since it exhibits at the 
same time phase coherence and density wave orders. In 
the following we will denote this phase as the DW-SF 
phase. Moreover it is easy to show from equations Q 
and pT7)) that 2Ps = ^dw as is numerically checked in Ta- 
ble |T1 As the order parameter is here of dimension three 
and the lattice is of dimension two, we do not expect any 
transition to an ordered phase at finite temperaturci^ 

Paiva et alX have studied the ground state of FRHM 
on a honeycomb lattice a few years ago. They found 
a QPT from an anti-ferromagnetic phase at large cou- 
pling to a metallic phase at low coupling, the critical cou- 
pling strength being bounded by 4 < Uc/t < 5. We use 



finite-size scaling and larger system sizes L to improve 
the numerical accuracy and narrow down the region of 
this QPT. Spin wave theory applied to Heisenberg mod- 
els implies that the structure and pair structure factors 
at T = scale with the number of lattice sites N = 2L^ 
hkei^i^i^ 

2P, (N) = S'dw {N)^aN + bVN + c 

where a,b,c are [/-dependent nonnegative constants. In 
the disordered phase Sdvi{N) is expected to reach a con- 
stant finite value as N goes to infinity, meaning that the 
coefficients a and b should then vanish. In the ordered 
phase, a should be strictly positive so that both Ps and 
5dw diverge linearly with N signaling the emergence of 
density and phase coherence orders. Using system sizes 
as large as L = 15, and using the vanishing of coefficient 
a to define the onset for the DW-SF phase, we have been 
able to infer the critical interaction strength Uc to be in 
the range 5.0 < Uc/t < 5.1 (Fig. [3]). 

In the study by Paiva et al., the metallic phase ap- 
pearing at low U was not studied in detail. In particular 
the question of the metallic or semi-metallic nature of 
the system was not addressed. Calculating the spectral 
function A{lli) for different values of U (Fig. [3]), wc find 
that the system is always a semi-metal when it is not in 
an ordered phase. The density of states drops around 
the Fermi level (located at w = 0) for [//t < 5 but with- 
out forming a gap. On the contrary, we observe a tiny 
metallic peak at the Fermi level. This peak is a finite-size 
effect due to the four states per spin located exactly at 
the Fermi level (in the non-interacting limit) when the 
system size is a multiple of three. On the contrary, us- 
ing sizes that are not multiples of three, we do observe 
a small gap. Both this gap and the peak are finite-size 
effects that are reduced when we increase the size of the 
system. We then conclude that A{uj) is zero (or very 
small) only at the Fermi level but without the formation 
of a gap. This is the signature of a semi-metallic phase. 
Indeed, a metal would be signaled by a persistent peak 
at the Fermi level (or at least a large non-zero density). 
The transition to the DW-SF ordered phase is signaled 
by the opening of the gap in A{uj) for U/t > 5, which 
corresponds to the value for the transition previously ob- 



n/t 


P 


S'dw/2 






0.9202 


1.0 
1.5 


1.125 ± 0.005 
0.3356 ± 0.0004 


1.127 ± 0.001 
10.5 ± 0.1 



TABLE I: Comparison of Ps and 5'dw/2 for L = 12, pt = 20, 
U/t — 3, and different values of /i/t. At half-filling, those 
quantities are equal within statistical error bars as a conse- 
quence of the SU(2) pseudo-spin symmetry of the FAHM. 
S'dw and Ps are small because U < Uc and the system is in 
its semi-metallic phase. This symmetry is broken when ^ 7^ 
and this is confirmed by the numerical data showing that the 
two quantities are indeed unequal. Sdw remains small but Ps 
is large due to the presence of quasi-long-range order. 
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FIG. 3: (color online) Scaling of the density wave structure 
factor Sdw with lattice size L at half-filling (the total number 
of lattice sites is A'^ = 2L^). The dashed lines are a fit of 
the form Sd„/N — a + b/^l4 + c/N . Close to or above the 
transition [U/t^ 5.0), the coefficients a and h take on finite 
positive values implying that both density and phase coher- 
ence orders emerge in the thermodynamic limit N —> <x. As 
it is seen, Sdw/N then essentially scales linearly with 1/V7^ 
and achieves the finite value a when N —> oo. Below the tran- 
sition {U /t ^ 5), the coefficients a and 6 vanish, meaning that 
the system reaches its disordered phase in the thermodynamic 
limit N ^ oo. As it is seen, Sdv,/N then essentially scales as 
and goes to zero when A'' oo. The QPT point is thus 
signaled by the vanishing of the coefficient a, from which we 
can infer that the critical interaction strength lie in the range 
5.0 < Uc/t < 5.1. 



tained by the finite-size scaling analysis of Sdv 



III. DOPING AWAY FROM HALF-FILLING 

At zero temperature, when the FAHM is doped away 
from the DW-SF ordered phase obtained at half-filling 
when U > Uc, say by increasing p from 1, we expect 
the density order to disappear and the phase coherence 
order to persist. However, one also expects phase coher- 
ence to establish throughout the sample when the system 
is doped away from the semi-metallic phase obtained at 
half-filling when U < Uc- Indeed in this case the Fermi 
surface is no longer limited to isolated points and BCS 
pairing becomes possible. Therefore, we expect the phase 
coherence order to establish at zero temperature for all 
values of the interaction U as soon as p 7^ 1. With an 
order parameter of dimension two (a phase gradient pic- 
tured as a vector lying in the xy-plane), the system un- 
dergoes a Berezinskii-Kosterlitz-Thouless (BKT)^i^i2^ 
transition at some critical temperature Tc, leading to a 
quasi-long-range phase order, i.e. a superfluid phase, at 
T < Tc before the appearance of the Bose-Einstein con- 
densate at T = 0. 



'>oU/t = 




FIG. 4: (Color online) Spectral function A{u) at half-filling 
(p = 1) for different values of the interaction strength U . The 
lattice size is L = 9 and fit = 10. The Fermi level is located at 
uj = G. For U/t < 5, the system is a semi-metal as witnessed 
by the dip around the Fermi level. The non- vanishing density 
of states at the Fermi level is due to finite-size effects (see 
Fig. [21). For U/t > 5, a. gap opens as the system enters the 
DW-SF ordered phase. The small peaks situated at \iu\ ~ 2.5 1 
are also a result of finite-size effects. 



According to mean-field theory^, a superconductor 
exists anywhere away from half-filling, albeit the su- 
perconducting gap function or, equivalently, (A|^), de- 
cays exponentially with respect to l/(Uy^p — 1) in the 
BCS regime. In their previous studyi^, Su et al. com- 
pared DQMC results to RPA calculations and showed 
that there is a so-called BCS-BEC crossover extend- 
ing from small to large values of the interaction when 
the system is off half-filling. When U is increased, the 
ground state of the system evolves continuously from a 
BCS state (where fermions with opposite spins form loose 
pairs of plane waves with opposite momenta) to a BEC 
of bosonic molecules (where fermions with opposite spin 
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FIG. 5: (Color online) Evolution of the pair structure factor 
Ps as a function of the inverse temperature f3t for several lat- 
tice sizes L. The interaction strength has been fixed a,tU = 3t 
and the total average fermionic density at p = 1.1. The 
dashed lines are fits using the 3-parameter function F{j3t), 
eq. (|ll|l . A plateau is reached when pt is much greater 
than the energy gap induced by finite-size effects between the 
ground state and the first excited state. As can be seen, the 
plateau is reached at larger pt when the lattice size increases. 
It is also reached at larger fit when p — > 1 (not shown) . 



form tightly-bound pairs) . We have extended their study 
to larger lattices (up to L = 15) and lower temperatures 
(up to f3t = 20) and we have also analyzed new obscrv- 
ables. 

We first studied the behavior of the pair and density 
wave structure factors, Ps and ^dw, away from half- filling. 
To do this, wc first need to obtain the low-tcmpcrature 
limit of these quantities by decreasing the temperature 
until we observe a plateau signaling that wc have reached 
the T = limit (Fig. [5]). To extract the plateau value, 
we have used the 3-paramcter function 



(11) 



to fit our numerical data Ps{f3t). The plateau value 
lim^^oo Ps is then approximated by u. We have also ob- 
served in our numerical simulations that this plateau is 
reached at lower and lower temperatures as we approach 
half-filling. This is because the BKT critical temperature 
Tc goes to zero like 1/| ln5p\ a,s 6p ~ \1 — p\ 0^^ and 
lower temperatures are required to achieve phase coher- 
ence. 

Fig. [6] shows how Pg and Sdw scale with the number 
of lattice sites N. For each chosen lattice size L and 
fermionic density p, we have run our simulations for the 
lowest temperature that could be numerically achieved. 
The temperature range that we have been able to explore 
was up to pt = 20. As expected always goes to 
zero and Pg always extrapolates to a non-zero value. We 
can then conclude, from direct measurement, that the 
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FIG. 6: (Color online) Evolution of the pair and density wave 
structure factors Ps and function of the number of 

lattice sites A'^ for different total average fermionic densities 
p. The interaction strength has been fixed at U = 3t. Full 
symbols have been obtained for inverse temperatures up to 
/3t — 20 (see text). Open symbols for Ps are the plateau 
values at T = as extracted from the fits in FigO The 
density wave structure factors Sd-w always go to zero as the 
system size L = \J N/2 tends to infinity whereas the phase 
coherence ordering signal Ps never vanishes. The dashed lines 
are guides to the eyes. For the same parameters at half-filling 
the system would be semi-metallic and Sdw and Ps would both 
vanish. 



BEC at zero temperature always appears as soon as the 
system is doped away from half-filling. Even with the 
smallest doping that we have been studying {p = 1.05, 
5% doping), we have observed a clear persistence of the 
phase coherence ordering in the large size limit. 

To observe the molecule formation along the BCS-BEC 
crossover, we have studied the density of on-site pairs 



(12) 



In the non-interacting limit {U/t 0), spin- up and 
spin-down particles are uncorrelated. Hence (n^in^) = 
(nil) (nil) = Pi Pi- Since we consider here equal spin 
populations p^ = pi ^ p/2, we find pp = p^. In the 
molecular limit {U/t oo), fermions can only exist in 
pair at a site. Hence {ni^riii) = {rii^) = p^ and p-p ~ p^. 
In see [71 we have plotted the rescaled density of on-site 
pairs: 



Pp = 



Pp - Pt 
P\- P\ 



(13) 



as a function of U/t. The BCS-BEC crossover is nicely 
evidenced by the smooth evolution of this rescaled quan- 
tity between the two limits /5p = and /jp = 1 as the 
interaction is increased. 
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FIG. 7: (Color online) Evolution of the rescaled density pp 
of on-site pairs, eq. (|13|l . as a function of the interaction 
strength U/t for two different total average fermionic den- 
sities p. The system size has been fixed at L = 9 and the 
inverse temperature is j3t = 10. In the non-interacting limit 
{U/t 0), spin-up and spin-down particles are uncorrelated, 
hence {nifUii) — {nii){nii) = PiPj. = P| for equal spin 
populations. In this case Pp = 0. In the molecular limit 
{U/t — » oo), fermions can only exist in pair at a site, hence 
(riiinii) = (nil) = p|. In this case pp = 1. The BCS-BEC 
crossover is characterized by the smooth evolution of pp be- 
tween these two limits and 1 as the interaction strength is 
increased. 



The second evidence for molecule formation along the 
BEC-BCS crossover comes from the evolution of the 
spectral function A{uj) when the interaction strength 
U (Fig. [5]) and the temperature T (Fig. [5]) are varied. 
At large interactions {U > 4), a clear gap is found at 
the Fermi level oj ~ provided the temperature is low 
enough, showing the formation of molecules. On the con- 
trary, when the interaction is weaker {U < 3), the gap 
does not open within the same range of temperatures. 
However, we observe that the value of A{uj) at the Fermi 
level u) = decreases when the temperature is lowered 
(Fig. [S]). We interpret this behavior as the precursor to 
the formation of a small BCS gap at very low temper- 
atures. This dip in A{uj) at the Fermi level is different 
from the one due to the vanishing of the non-interacting 
density of states at the Dirac points that was observed at 
half-filling in the semi-metal case. The Dirac dip is still 
present in the [/ < 3 cases for w < (Fig. [8]), showing 
that interaction strength is not large enough to strongly 
modify the structure of the Fermi sea, except very close to 
the Fermi level. This is characteristic of the BCS case. 
On the other hand, the Dirac dip disappears at strong 
interactions (Fig. [51 bottom), showing now that the orig- 
inal Fermi sea structure has been completely modified by 
interactions. 

A nice feature of the strongly-interacting regime is 
the existence of two very different energy scales. One 




(a) 





FIG. 8: (Color online) Evolution of the spectral function A{uj) 
as a function of the interaction strength U at density p = 1.2, 
inverse temperature j3t = 12 and lattice size L = 9. When 
(7 = 0, the chemical potential is numerically found to be 
p/t = 0.8768, locating the Dirac points in the residual gap 
(due to finite-size effects and temperature rounding) around 
uj/t = —1. The fact that the density of states vanishes linearly 
with u around u/t = —1 also supports this identification of 
the location of the Dirac points. As U is increased, a dip 
develops in the spectral function at the Fermi level (located 
at = 0) and the BCS-BEC gap eventually opens while the 
Dirac points are gradually destroyed. 



corresponds to the formation of tightly-bound pairs 
(molecules) and is typically of the order of U itself. The 
second corresponds to the emergence of phase coherence 
between these pairs and is of the order of the hopping 
parameter for pairs, typically t} /U These two energy 
scales are clearly identified by comparing the evolution 
of Ps and /Op when the temperature is varied, see Fig. [TOl 
We thus can conclude that, even at \J /t = 3, we ob- 
serve the formation of pairs before the emergence of phase 
coherence. To investigate this phenomenon further, we 
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15 



FIG. 9: (color online) Evolution of A{uj) as a function of 
inverse temperature pt at p — 1.2, interaction strength U = 2t 
and lattice size L — 9. As the temperature is lowered, a dip 
develops in the spectral function at the Fermi level located at 
uj — 0. Eventually a gap opens when the temperature is low 
enough (not shown). The gap opening at the Fermi level is 
obtained even at weak interactions, a situation characteristic 
of the existence of a small BCS gap. 



show in Fig. [TT] the pair Green's function (|3]) as a func- 
tion of distance for different temperatures. There is a 
range of temperatures (0.1 < pt < 5) where the pair 
Green's function is clearly decreasing exponentially with 
distance (up to some boundary effects). This means that 
no phase coherence is achieved and the system is in a 
disordered regime. In other words, the corresponding 
temperatures are above the BKT transition temperature 
Tc- For this same temperature range, pp has already 
reached its zero-temperature limit (Fig. [TU)) . This is a 
clear evidence for the existence of preformed pairs which 
will eventually develop quasi-long-range phase coherence 
at a much lower temperature. For temperatures T < T^, 
the Green's function should decay algebraically with dis- 
tance with an exponent rj = T/{'iTc)i^ For l3t > 10, the 
pair Green's function behavior is consistent with a power- 
law decay, but it is difficult to extract the corresponding 
exponent due to finite-size effects. 



IV. CONCLUSION 

We have studied the Hubbard model on a honeycomb 
lattice with attractive interactions. At half-filling, build- 
ing up upon previous existing studies, we have used the 
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FIG. 10: (Color online) Evolution of the pair structure fac- 
tor Ps (circles) and the rescaled density of on-site pairs pp 
(squares) as a function of the inverse temperature /3t at inter- 
action strength U = 3t. The total average fermionic density 
is set at p = 1.5 and the system size is L = 12. Two dif- 
ferent energy scales are clearly identified as Ps, signaling the 
emergence of phase coherence, saturates at /3t ~ (7 /t whereas 
Pp, signaling the molecule formation, saturates at j3t « t/U. 
We recover here (in dimensionless units) the two energy scales 
/U and U, typical of the emergence of phase coherence and 
of the formation of tightly-bound pairs. 



mapping onto the FRHM to show that there is a quan- 
tum phase transition at T = between a disordered 
phase and a DW-SF phase exhibiting crystalline as well 
as superfluid orders. The critical interaction strength at 
which this QPT takes place is accurately bounded by 
5.0 < Uc/t < 5.1. We have also shown that, before the 
transition, the system is semi-metallic and that the inter- 
actions do not markedly change the nature of this phase. 
Away from half-filling, within our numerical accuracy, 
the system seems to become superfluid, even for arbi- 
trary small values of the doping. We have elucidated the 
presence of the BCS-BEC crossover by looking at several 
quantities, especially the one-particle density of states. 
We have clearly evidenced, for strong enough interac- 
tions, the existence of two different energy scales, one for 
the formation of the pairs and one for the emergence of 
phase coherence (the BKT transition), which is typical 
of the strongly interacting regime. 

For weak interactions, both at and away from half- 
filling, we have observed that the spectral function A{llj) 
is qualitatively the same as in the non-interacting case. 
Only the states close to the Fermi level are affected by 
those weak interactions. As there are no available states 
in the half-filled case close to the Fermi level, the interac- 
tions hardly play a role and the system remains a semi- 
metal (at half-filling) up to f/ = 5t. It is only when the 
interactions are strong enough to destabilize the Fermi 
sea and form tightly-bound pairs that the system enters 
a different phase. In this case, the description in terms 
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FIG. 11: (Color online) Evolution of the pair Green's function 
as a function of distance for different temperatures. The to- 
tal average fermionic density is set at p — 1.5, the interaction 
strength at U = 3t and the lattice size is L = 12. The vertical 
axes are plotted in logarithmic scale while the horizontal axes 
are plotted with linear (top) and logarithmic (bottom) scales. 
For large site separation \i — j\, we observe a transition from 
an exponential decay (linear behavior in the log-linear plot) at 
high temperature to a weak algebraic decay (linear behavior 
in the log-log plot) at low temperature. This is the signature 
of the BKT transition where the system leaves the disordered 
phase to enter a phase with quasi-long-range order as the tem- 
perature is lowered. However, due to limited system size, the 
weak algebraic decay of the pair Green's function is difficult 
to infer unambiguously. 



of individual fermions and plane- wave states is no longer 
relevant. 



We further observe that the BCS and the semi-metal 
regimes are two phases sharing some common features. 
Indeed, in both phases, interactions are not strong 
enough to substantially modify the Fermi sea structure 
except around the Fermi level. This is reflected in the 
fact that the Dirac dip in A{uj) is always clearly visi- 
ble in these cases. By the same token, the molecular 
superfluid phase (EEC) and the DW-SF have in com- 
mon that the description in term of individual fermions 
is meaningless. Indeed, for both phases, the fermionic 
excitations are gaped and the Dirac dip in A{oj) has dis- 
appeared. Close to half-filling, we then observe the BCS- 
BEC crossover to happen for interaction strengths close 
to the value of the QPT at half-filling, i.e. U = 5t. 
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